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1  Introduction 


The  primary  objective  of  this  project  has  been  to  develop  effective  ultrasonic 
tomographic  algorithms  that  are  applicable  to  the  endoluminal  geometry 
found  in  prostate  imaging.  The  method  of  reconstruction  is  that  of  inverse 
scattering  in  which  an  unknown  scatterer  is  recovered  from  the  knowledge 
of  the  fields  scattered  by  the  object  in  an  appropriately  designed  scattering 
geometry.  It  was  mentioned  in  our  previous  reports,  that  inverse  scattering 
is  capable  of  overcoming  the  two  major  limitations  of  conventional  ultra¬ 
sound  imaging,  namely,  its  non-quantitativeness,  subjectivity  and  sensitivity 
to  speckle  noise  arising  from  multiple  scattering  of  acoustic  wave  radiation. 
However,  numerous  challenges  remain  before  inverse  scattering  ultrasonic 
tomography  becomes  feasible  in  real-world  applications.  It  is  well-known 
that  ultrasound  tomography  in  inverse  scattering  is  mathematically  com¬ 
plex,  computationally  intensive,  and  like  any  inverse  solution  of  a  physical 
problem,  ill-posed.  [1,  3]  In  the  previous  report,  we  demonstrated  how  the  ill- 
posedness  could  be  overcome  by  using  regularization,  appropriate  data  sets 
and  reconstruction  techniques.  Toward  that  a  multi-wave,  multi-frequency 
data  set  was  used  in  Tikhonov  regularization  with  a  particularly  chosen  L-2- 
norrn  in  the  object  error.  Moreover,  the  reconstructions  were  performed 
by  stepping  up  in  grid  sizes  instead  of  computing  with  the  full  grid  size 
all  at  once.  We  demonstrated  improved  convergence  using  frequency  hop¬ 
ping  multigrid  inversion,  and  showed  that  it  was  possible  to  overcome  the 
problem  of  the  local  minima  and  improve  inversion  performance  using  this 
approach.  Furthermore,  the  regularizing  character  inherent  in  multi-gridding 
was  also  demonstrated  in  the  previous  report.  A  further  significant  gain  in 
computational  speed  was  achieved  by  using  the  adjoint  fields  to  calculate  the 
gradient  of  the  objective  function  eliminating  thereby  the  computationally 
intractable  task  of  having  to  determine  the  full  Jacobian  matrix.  The  ap¬ 
plication  of  inverse  scattering  ultrasonic  tomography  in  quantitative  tissue 
imaging  was  convincingly  demonstrated  in  our  last  progress  report  via  the 
reconstructions  2-D  and  3-D  Shepp- Logan  tissue  phantoms  by  nonlinear  con¬ 
jugate  gradient  method  on  the  basis  of  multi- incident,  multi-frequency  data 
and  multi-grid  reconstruction.  In  addition,  the  feasibility  of  the  technique 
was  further  tested  by  reconstructing,  in  addition  to  the  tissue  phantoms,  a 
stylized  endoluminal  geometry,  and  its  was  done  in  two  different  scattering  ge¬ 
ometries  -  the  usual  exterior /exterior  geometry  and  a  novel  exterior/interior 
or  scmi-endoluminal  scattering  configuration  which  was  specially  designed 
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with  a  view  to  transurethral  application. 

In  the  final  year  of  this  project,  we  have  continued  to  address  the  is¬ 
sue  of  accelerating  inversion  performance.  This  work  has  focused  on  imple¬ 
mentation  of  a  preconditioning  scheme  for  accelerating  convergence  of  the 
conjugate  gradient  minimization  algorithm,  developing  and  coding  a  fast 
contrast-source  inversion  algorithm  to  provide  an  improved  initial  guess  at 
the  scattering  object,  and  design  and  implementation  of  a  scheme  for  par¬ 
allel  computation  of  the  forward  model  using  the  Message  Passing  Interface 
(MPI).  We  have  also  recently  developed  a  collaboration  with  a  company 
producing  a  commercial  scanner  for  ultrasound  tomography  (in  the  exte¬ 
rior/exterior  geometry)  and  have  begun  to  test  our  reconstruction  algorithms 
using  test  data.  We  have  also  developed  a  theoretical  analysis  of  noise  prop¬ 
agation  in  inverse  scattering  (attached  in  Appendix  A),  analyzed  acoustical 
wave  propagation  in  circular  ducts  of  various  wall  conditions  in  the  presence 
of  flow  (attached  in  Appendix  B),  and  investigated  transmitted  acoustic  field 
patterns  for  catheter-based  I  YUS  transmitters. 


2  Preconditioning 

As  regarding  the  solution  to  the  problem  of  the  computational  complexity  of 
inverse  scattering  tomography  in  imaging  the  prostate  (or  any  tissue  complex 
and/or  endoluminal  anatomy,  for  that  matter),  two  recourses  are  available 
in  general:  parallel  computing  and  designing  a  suitable  preconditioner.  Both 
are  discussed  in  this  report,  beginning  with  the  preconditioner  first.  However, 
in  view  of  the  complexity  of  the  matrix  mathematics  that  is  involved  in  the 
preconditioner  calculation,  only  the  basic  outlines  and  the  final  results  are 
given. 

The  preconditioner  described  here  is  based  upon  the  seminal  work  of 
Hohage.[7]  The  basic  idea  behind  the  preconditioner  is  to  combine  the  well- 
known  Lanczos  algorithm^ ]  for  calculating  the  eigenvalues  of  a  matrix  (via 
tri-diagonalization) ,  and  the  orthogonal  sequences  of  the  conjugate  directions 
that  are  generated  during  the  iterations  (NLCG  iteration  in  this  project). 
The  first  step  in  the  design  of  the  preconditioner  is  to  transform  the  varia¬ 
tional  problem  (detailed  in  the  previous  two  reports)  into  normal  equations. 
It  involves  transforming  the  objective  functional  (that  was  minimized)  into 
the  form:  =  F'n{l*)T F'n(l*)  +  ocnI.  F'  is  the  Frechet  derivative  of  the 

functional  with  respect  to  the  unknown  object  7,  and  an  is  an  iteratively 


3 


updated  regularization  parameter.  I  is  the  identity  matrix,  the  superscript 
T  denotes  the  transpose  operation,  and  n  is  the  stage  of  iteration. 

The  preconditioner  is  obtained  from  the  eigenvalues  and  eigenfunctions 
of  the  normal,  symmetrical  matrix,  F^*)T F'n{^*)n.  Let  {vj,\ j,  >  A2  > 
A3  •  •  •  }  be  the  eigenvalues  and  eigenfunctions.  Let  r  be  the  orthonormal- 
ized  conjugate  direction  vector  that  is  generated  within  the  CGNE  rou¬ 
tine,  and  q  the  vector  which  appears  in  the  object  (in  NLCG)  or  in  the 
Newton  step  (in  the  case  of  the  Newton  method)  update.  In  other  words, 
Xn+i  =  Xn  +  n  being  the  stage  of  iterartion.  The  Lanczos  tri¬ 

diagonal  matrix  is  built  from  r  and  q  vectors.  The  actual  computation  of 
the  spectrum  of  this  tri-diagonal  matrix  can  be  accomplished  by  any  of  the 
available  robust  method  for  the  purpose  such  as  the  QR-decomposition.  [4] 
The  preconditioner,  P,  takes  the  form: 


k 

Pnx  =  anx  +  ^  A j(x,  Vj)vj, 
3= 1 


from  which  it  follows  that 


3  Contrast-Source  Inversion  (CSI) 

In  our  last  report,  mention  was  made  of  the  CSI  method.  We  have  finished 
writing  the  MATLAB  code  for  multidimensional  CSI  reconstruction  using 
multifrequency,  multi-wave  data.  There  are  two  components  in  CSI,  namely, 
the  source  and  the  object.  Although  the  two  are  interdependent,  their  up¬ 
dates  can  be  performed  separately.  With  the  following  definitions  : 

scat  =  the  data  vector. 

C\  =  norm  (scat)2. 

7O  =  the  initial  estimate  of  the  object  obtained  by  weighted  backpropagating 
of  the  measured  data. 

7  =  the  reconstructed  object. 

Wjn  =  the  equivalent  source  in  the  j-th  incidence  and  in  the  n-th  stage  of 
iteration. 

Ujn  =  the  total  held  corresponding  to  Wjn. 


4 


u™c  =  the  incident  field. 

Pjn  =  the  data  discrepancy  in  the  j-th  incidence  and  in  the  n-th  stage  of 
iteration. 

r jn  =  the  object  discrepancy  in  the  j-th  incidence  and  in  the  n-th  stage  of 
iteration. 

Go  —  Green’s  function  in  the  computational  domain. 

Gs  =Green’s  function  in  the  detector  domain. 
gfn—  gradient  in  the  w-update. 

C  =  steplength  in  the  update  of  the  Polak-Ribiere  conjugate  directions  for 
w. 

Vjn  =  conjugate  directions  in  the  linear  conjugate  gradient  iteration  for  w. 
awn  =steplength  in  the  update  of  w. 

(yb  =  gradient  in  the  object-update. 

£yn=  steplength  in  the  update  of  the  Polak-Ribiere  conjugate  directions  for 
the  object. 

=  conjugate  directions  in  the  linear  conjugate  gradient  iteration  for  the 
object. 

ain  =  steplength  in  the  update  of  the  object. 

*  =  the  operation  of  convolution. 

CD  =  the  computational  domain. 

DD  =  the  detector  domain. 

(-,  -)y  =  the  inner  product  on  Y. 

Y2’  =sum  over  j  and  n. 
z—  the  complex  conjugate  of  z. 

our  implementation  of  CSI  is  as  follows: 

The  Source  Update: 
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C^'inn. 


Cr2,n- 1  —  It? 


-1« 


.mc|2 
3 


'U'j,n+ 1  =  Uj  ~\~  G D  *  Wjn. 

1  Wjn  H“  QjnVjn* 
^jn  Sjfn  “1“  ^wn^wn- 

=  Re  E'  <  9jn>  9jn  ~  9j,n- 1  >CT 
<9f,n-l,9j,n-l>CD 

<  afn,  Vjn  >CD 


C\ (norm(Gsvjn)2  +  C2,n-i(norm(vjn  -  jn-i GDvjn)2' 
din  =  [Cl  G  sPj,n—l  +  —  GoJu-ijr^n-i]. 


The  Contrast  Update: 


9n  —  C2,n-1 


E,'  WjnU 


jn 


7n-l 


-  7'n 


Ej  lwjnh 

d)l  (j  T  £ryndn_  1. 
ReJ2<  9l,9l-9l-i  >CD 
<  9n-l’9n- 1  >CD 

7??+l  T?i  T  OLlndn. 


The  Polak-Ribiere  steplength,  cc7n,  is  a  highly  involved  expression  contain¬ 
ing  the  current  updates  in  both  the  source  and  the  contrast  quantities,  and 
its  explicit  expression  is  not  given  here.  We  applied  the  CSI  algorithm  to 
multifrequency  data,  but  upsampling  in  frequency. 

The  backpropagated  CSI  image  was  incorporated  into  the  NLCG  routine 
showing  a  slight  reduction  in  the  number  of  iterations.  This  suggests  that  a 
correctly  reconstructed  CSI  image  (which  is,  as  already  mentioned,  is  quite 
fast)  as  an  input  to  the  NLCG  iteration  would  lead  to  faster  convergence. 
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4  Acquiring  Experimental  Data 

Only  recently,  we  have  succeeded  in  gaining  access  to  a  commercial  scanner 
which  can  be  used  to  generate  experimental  data.  The  scanner  is  owned 
by  a  locally  based  breast  imaging  company,  Techniscan  Medical  Systems, 
Inc.  As  of  now,  Techniscan’s  scanner  operates  in  the  transmission  mode. 
However,  it  is  soon  to  be  modified  for  reflection  mode  imaging  which  is  what 
is  required  for  prostate  imaging.  We  have  already  acquired  some  trial  run, 
transmission  mode  data.  The  details  of  the  Techniscan  scanner  in  which 
the  data  were  acquired  is  shown  in  Figure  1.  [9] .  It  is  essentially  a  confocal 
geometry  in  which  two-dimensional,  slice-by-slice,  180°  data  can  be  collected 
in  the  temporal  as  well  as  in  the  frequency  domain.  The  data  that  we  acquired 
were  the  frequency  domain  data.  Figure  1.  Techniscan  Medical’s  ultrasonic 
scanner  geometry,  (a)  The  128  mm  x  16  mm  receiver  array.  There  are  6 
rows  and  160  columns,  i.e.,  960  elements  total.  The  central  two  rows  are  4.6 
mm  each.  The  next  outer  two  rows  are  1.9  mm  each,  whereas  the  rest  two 
rows  are  1.45  mm  each,  (b)  The  plan  view  of  the  scanner.  The  entire  system 
is  immersed  in  a  water  bath.  The  transmitter-receiver  system  rotate  together 
around  an  axis.  The  measurements  are  made  at  every  2°  of  roration.  (c)  The 
focussing  arrangement.  Both  the  transmitter  and  the  receiver  are  focussed 
with  a  focal  length  of  60  mm. 

Figure  1(a)  shows  the  details  of  the  scanner’s  receiver  array.  The  array 
has  960  receiving  elements.  The  length  of  the  array  is  128  mm  and  its  width  is 
16  mm  containing  6  rows.  The  central  two  rows  are  each  of  width  4.6  mm,  the 
adjacent  two  are  1.9  mm  wide  each,  whereas  the  outermost  2  rows  are  each 
1.45  mm  wide.  The  transmitter  which  has  the  length  of  140  mm  and  width  of 
16  mm,  is  made  from  a  single  piece.  Figure  1(b)  is  the  plan  view  of  the  scanner 
showing  the  organ  (scatterer)  between  the  transmitter  and  the  receiver,  the 
entire  arrangement  being  immersed  in  a  water  bath.  In  order  to  acquire  two- 
dimensional  slices  at  discrete  “levels” ,  that  is,  at  various  heights  in  the  object 
(a  prostate  phantom,  for  example),  the  transmitter  and  the  receiver  are  both 
focussed  (focal  length  =  60  mm),  as  shown  in  Figure  1(c),  which  displays  the 
view  of  the  elevation.  The  data  are  first  acquired  by  rotating  the  transmitetr- 
receiver  system  around  the  scatterer.  The  measurements  are  made  at  every 
2°  of  rotation.  Because  of  the  focusing  arrangement,  approximate  2-D  data 
at  any  pre-set  level  can  be  obtained  by  simply  summing  over  the  six  bins  of 
the  receiver  array  in  Figure  1(a),  that  is,  along  the  width  of  the  array.  The 
time  domain  data  are  acquired  by  chirping  the  transmitter.  These  are  then 
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Fourier  transformed  into  the  frequency  domain  data.  The  clock  rate  is  33 
MHz,  and  the  FFT  is  performed  after  1536  time  intervals. 

5  Parallelization  of  the  Forward  Model 

For  realistic  imaging  geometries,  it  is  expected  that  anywhere  from  64-256 
incident  waves  and  a  corresponding  number  of  detectors  will  be  used  in  the 
image  reconstruction  system. 

6  Noise  Propagation  in  Inverse  Scattering 

The  propagation  of  noise  from  the  detectors  to  the  image  via  the  recon¬ 
struction  algorithm  was  investigated.  The  algorithm  in  this  project  was: 
V7$  =  Re  —  0.  Expressions  for  the  covariance  matrix  was  derived. 

The  covariance  reduced  to  that  of  the  X-ray  CT  under  the  assumptions  of 
linear  operator  and  real  data.  [5]  The  analytical  covariance  matrix  was  further 
specialized  to  the  limit  of  the  Born  approximation.  The  variance  in  the  Born 
image  pixels  was  found  to  be  constant,  while  the  covariance  at  a  pixel  was 
highly  peaked  and  decayed  oscillatorily  away  form  the  pixel.  These  Endings 
were  in  agreement  with  the  results  obtained  by  the  direct  Born  calculations. 
Further  details  are  presented  in  ??. 

7  Acoustic  Emission  from  an  IVUS  Trans¬ 
ducer 

In  the  last  progress  report,  some  calculations  were  presented  regarding  the 
characteristics  of  the  fields  emitted  by  a  catheter-based  miniature  transmit¬ 
ter.  For  transurethral  imaging  of  the  prostate  gland,  the  knowledge  of  the 
characteristics  of  such  transmitters  is  of  interest.  Some  results  of  numerical 
computations  are  shown  in  Figures  2.  In  this  figure,  a  miniaturized  emitter 
is  situated  on  the  cylindrical  surface  of  the  catheter,  the  diameter  of  which 
is  3  F.  The  emitter  has  a  total  angular  width  of  0.375°.  The  wavenumber 
k0  =  50.  For  the  Robin  or  the  impedance  catheter,  the  wall  admittance  in 
the  plots  shown  is  real  and  its  value  is  15.  The  plots  are  representatives  of 
the  emitted  fields. 


Figure  2(a)  results  if  the  catheter  surface  is  a  Neumann  surface,  whereas 
Figure  2(b)  shows  the  results  for  a  Robin  or  impedance  catheter.  The  fields 
emitted  by  a  transmitter  in  an  unbounded  space  are  shown  in  Figure  2(c). 
A  comparison  between  Figures  2(a)  and  2(b)  with  Figure  2(c)  clearly  shows 
that  the  fields  of  the  catheter-based  transmitters  are  more  directional  than 
the  transmitter  in  the  unbounded  medium.  The  latter  is  an  isotropic  emit¬ 
ter  as  opposed  to  the  catheter-based  ones.  The  objective  is  to  understand 
if  a  two-dimensional,  slice-by-slice  imaging  is  possible  with  the  miniature 
transmitters.  Further  studies  are,  however,  needed  for  a  definitive  answer. 

8  Key  Research  Accomplishments 

The  central  objective  of  the  project  was  to  develop  methods  for  imaging  in¬ 
ternal  anatomies  with  an  emphasis  on  imaging  tissue  complexes  such  as  the 
prostate  gland  by  quantitative  ultrasound  so  that  some  of  the  deficiencies 
of  conventional  ultrasound  imaging  could  be  removed.  Inverse  scattering  in 
ultrasonic  frequencies  was  proposed  as  the  method  of  choice.  Two  key  issues 
are  involved  in  the  implementation  of  any  ultrasonic  inverse  scattering  algo¬ 
rithm.  These  are:  (1)  efficient  solutions  of  the  Lippmann-Schwinger  integral 
equation  of  scattering  -  the  so-called  forward  problem,  and  (2)  the  speed  of 
the  numerical  inversions.  The  latter  in  turn  includes  the  problem  of  fast 
computation  of  the  gradient,  negotiation  of  the  local  minima,  and  general 
imaging  protocols  for  gaining  advantage  in  speed.  The  project  accomplish¬ 
ments  are  summarized  below. 

•  Fast  and  efficient  codes  for  the  solution  of  the  Lippmann-Schwinger 
equation  were  written  using  CG-FFT.  These  were  checked  against  ana¬ 
lytical  benchmark  solutions.  Our  numerical  results  were  found  to  have 
a  negligible  error  of  about  only  0.4  percent.  It  was  a  convincing  demon¬ 
stration  of  the  efficacy  and  accuracy  of  the  numerical  codes  which  were 
developed  in  this  project  for  solving  the  forward  problems  in  2-D  and 
in  3-D. 

•  For  inversions,  multi-frequency  data  were  used  in  order  to  minimize 
the  problem  of  the  local  minima.  Moreover,  the  adjoint  fields  which 
were  used  in  the  calculation  of  the  gradient  of  the  objective  functional, 
led  to  a  significant  reduction  in  computations.  Only  one  extra  forward 
problem  (per  incidence)  needed  to  be  solved  in  this  technique.  The 
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gradient,  calculated  by  the  method  of  adjoint  fields,  resulted  in  almost 
identical  results  when  compared  against  the  brute  force  finite  difference 
calculations. 

•  The  multigrid  bootstrapping  approach  yielded  relatively  speedier  re¬ 
constructions. 

•  An  exterior- interior  or  semi-endoluminal  configuration  was  introduced 
and  investigated  for  intra-organ  recosntructions  during  the  course  of 
this  project.  It  is  an  entirely  novel  scattering  geometry  which  is  ex¬ 
pected  to  be  useful  for  imaging  any  internal  tissue  anatomy  as  well  as 
imaging  the  lumen  of  a  vessel. 

•  The  feasibility  of  the  imaging  procedures  that  were  proposed  and  devel¬ 
oped,  including  the  new  semi-endoluminal  configuration,  were  convinc¬ 
ingly  demonstrated  by  numerically  reconstructing  Shepp-Logan  tissue 
phantoms  and  stylized  endoluminal  phantoms  in  both  2-D  and  in  3-D. 
In  addition,  the  strength  of  the  scatterer  as  high  as  7 max  =  0.3,  that 
is,  allowing  for  as  much  as  30  percent  variation  in  the  speed  of  sound 
relative  to  the  background.  It  is  well  beyond  the  Born  limit. 

•  As  an  important  contribution,  we  have  outlined  a  preconditioner  that 
was  first  proposed  by  Hohage  for  the  Newton  methods  using  normal 
equations.  The  preconditioner  led  to  orders  of  magnitude  gain  in  the 
speed  of  the  inversions.  We  have  proposed  the  use  of  Hohage’s  precon- 
clitioner  for  the  NLCG  iterations,  as  were  done  in  this  project. 

•  We  have  written  MATLAB  code  for  rapid  inversions  by  the  Contrast- 
Source-Inversion  (CSI)  technique  that  uses  the  perconditioned  gradient 
and  only  linear  CG.  All  expressions  in  CSI  are  analytical,  and  the 
method  yields  quite  rapid  inversions.  However,  our  purpose  is  to  use 
the  basic  CSI  in  a  bootstrapped  frequency  and  upsampling  to  accelerate 
convergence  and  improve  accuracy  of  our  reconstructions. 

•  An  arrangement  was  made  with  Techniscan  Medical  Systems,  Inc.,  Salt 
Lake  City,  Utah,  for  using  its  scanner  for  acquiring  experimental  patient 
data.  Techniscan  specializes  in  breast  tumor  imaging.  Both  temporal 
and  spectral  data  can  be  obtained.  We  have  already  acquired  frequency 
domain  data  acquired  in  vivo  in  the  human  breast  and  have  begun  tests 
of  our  algorithm  using  this  data. 
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•  A  scheme  of  a  parallel  computation  for  accelerating  the  solution  of 
the  inverse  problem  has  been  developed  and  implemented  using  MPI. 
This  approach  facilitates  the  use  of  widely  available  computational  clus¬ 
ters  (MIMD)  supporting  this  standard  protocol.  While  communication 
overhead  limits  performance  for  small  grid  sizes,  this  becomes  relatively 
unimportant 

•  Numerical  results  have  been  obtained  for  the  transmitted  acoustic  fields 
of  catheter-based  IVUS  transducers. 

•  Novel  theoretical  results  on  noise  propagation  in  inverse  scattering  us¬ 
ing  the  adjoint  held  gradient  have  been  developed. 

•  Peripheral  to  the  investigation  of  scmi-endoluminal  scattering,  the  char¬ 
acteristics  of  propagation  of  acoustic  waves  in  cylindrical  waveguides 
for  various  types  of  wall  conditions  were  investigated  both  theoretically 
and  numerically. 


9  Reportable  Outcomes 

The  accomplishments  cited  above  have  been  translated  into  an  invited  lec¬ 
ture,  a  refereed  journal  publication,  and  one  paper  submitted  to  a  peer- 
reviewed  journal,  and  one  paper  in  preparation  for  submission  to  a  peer- 
reviewed  journal: 

•  DN  Ghosh  Roy,  John  Roberts,  Matthias  Schabel  and  SJ  Norton,  Noise 
propagation  in  linear  and  nonlinear  inverse  scattering,  J.  Acous.  Soc. 
Am.,  vol.  125,  May  2007. 

•  DN  Ghosh  Roy,  John  Roberts  and  Matthias  Schabel,  Acoustic  wave 
propagation  in  a  cylindrical  waveguide  in  the  presence  of  how,  submit¬ 
ted  to  J.  Sound  and  Vibration. 

•  Matthias  Schabel  and  DN  Ghosh  Roy,  Two  and  three  dimensional  ul¬ 
trasound  tomography  in  exterior  and  semi-endoluminal  geometry,  in 
preparation  for  submission  to  Physics  in  Medicine  and  Biology. 

•  DN  Ghosh  Roy,  John  Roberts,  Matthias  Schabel  and  SJ  Norton,  The 
Method  of  Adjoint  Field  in  Inverse  Scattering  of  Plane  Waves ,  Invited 
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Talk  delivered  at  the  Acoustical  Society  of  America  Annual  Convention, 
Providence  RI,  June  2006. 


10  Conclusions  and  Future  Work 

This  project  set  an  ambitious  research  plan;  while  not  every  proposed  element 
of  the  work  detailed  in  the  initial  grant  submission  has  been  fully  achieved, 
we  have  made  substantial  progress  on  a  number  of  fronts.  In  particular, 
we  have  developed  and  implemented  code  for  performing  full- wave  acoustic 
inversion  that  is  functions  in  2D  and  3D  and  for  both  the  exterior-exterior 
and  exterior-interior  source/detector  geometries.  We  have  demonstrated  the 
feasibility  of  frequency  hopping  as  a  means  of  eliminating  the  need  for  em¬ 
pirical  regularization  of  the  inversion  algorithm  and  have  demonstrated  the 
ability  to  reconstruct  large  objects  in  the  presence  of  simulated  noise  with 
excellent  fidelity  (up  to  grid  sizes  of  256x256  in  2D  and  40x40x40  in  3D). 
Computation  time  for  reconstruction  remains  a  challenge. 
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Figure  1:  Figure  1.  The  geometry  of  the  ultrasound  scanner  for  acquiring  nl- 
trasonically  scattered  data  in  time  and  frequency,  (a)  The  receiver  array:  128 
mm  x  16  mm,  960  total  elements,  (b)  The  plan  view  of  the  scanner  showing 
the  scatterer  and  the  transmitting- receiving  arrangement.  The  entire  set  np 
is  immersed  in  a  water  bath.  The  transmitter  and  receiver  rotate  as  a  solid 
body  around  the  axis  shown.  Data  collection  is  at  every  2°  of  rotation,  (c) 
The  focusing  geometry.  Both  the  transmitter  and  the  receiver  are  focused 
with  a  focal  length  of  60  mm.  The  transmitter  has  a  lens  for  focusing  while 
the  receiver  is  surface  is  curved. 
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Angular  Distribution:  No  Cylinder  Angular  Distribution:  Neumann  Cylinder 

0o=O.  1  00/2  k=  1  5  0 0=0.375/2  k=50 


Figure  2:  Plots  showing  the  angular  variation  of  the  held  emitted  by  a 
catheter-based  miniature  emitter,  (a)  A  Neumann  catheter.  k0  =  50,  and 
0o  =  0.185°.  (b)  A  Robin  catheter,  a  =  15.  k0,  0Q  as  in  (a),  (c)  Emitter 
radiating  in  an  unbounded  space.  Other  parameters  same  as  in  (a). 
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The  propagation  of  noise  from  the  data  to  the  reconstructed  speed  of  sound  image  by  inverse 
scattering  within  the  framework  of  the  Lippmann-Schwinger  integral  equation  of  scattering  is 
discussed.  The  inversion  algorithm  that  was  used  consisted  in  minimizing  a  Tikhonov  functional  in 
the  unknown  speed  of  sound.  The  gradient  of  the  objective  functional  was  computed  by  the  method 
of  the  adjoint  fields.  An  analytical  expression  for  the  inverse  scattering  covariance  matrix  of  the 
image  noise  was  derived.  It  was  shown  that  the  covariance  matrix  in  the  linear  x-ray  computed 
tomography  is  a  special  case  of  the  inverse  scattering  matrix  derived  in  this  paper.  The  matrix  was 
also  analyzed  in  the  limit  of  the  linearized  Born  approximation,  and  the  results  were  found  to  be  in 
qualitative  agreement  with  those  recently  reported  in  the  literature  for  Born  inversion  using  filtered 
backpropagation  algorithm.  Finally,  the  applicability  of  the  analysis  reported  here  to  the  obstacle 
problem  and  the  physical  optics  approximation  was  discussed.  ©  2007  Acoustical  Society  of 
America.  [DOI:  10.1121/1.2713671] 
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I.  INTRODUCTION 

Inverse  scattering 1-4  is  an  intensely  researched  area  in 
current  applied  mathematics,  and  has  wide  applications  in 
almost  all  branches  of  modern  science  and  technology,  from 
medical  imaging  and  nondestructive  evaluation  of  materials 
to  problems  in  astronomy.  The  general  objective  is  to  recover 
an  unknown  physical  object  (scatterer  and/or  inhomogeneity) 
from  a  set  of  appropriately  defined  data.  The  data  in  inverse 
scattering  consist  of  fields  scattered  by  the  unknown  object 
under  wave  excitation.  The  framework  of  analysis  in  this 
work  is  the  famous  Lippmann-Schwinger  integral  equation  of 
scattering. 1-4  Let  T  denote  the  integral  operator  of  scattering. 
If  the  scatterer  is  “weak,”  then  the  uniform  norm  of  T, 
||7]|oo^l,  where  ||7]|O0  =  (l/2)(A:0a)2||  1 -n)^.  1  -n  is  the  devia¬ 
tion  of  the  scatterer  refractive  index  n  relative  to  that  of  the 
background  which  was  assumed  to  be  unity,  k()=2Tr/k  is  the 
wave  number  of  the  exciting  wave  which  is  of  wavelength  X. 
a  is  the  dimension  of  the  scatterer.  Moreover,  ||x(t)||oo 
=  max|x(r)|  over  a  certain  interval  of  t  and  x  in  a  suitable 
space  X.  In  this  case,  the  inverse  scattering  problem  can  be 
linearized  by  considering  only  the  single  scattering  events. 
This  is  the  well-known  Born  and  Rytov  approximation/  In 
these  approximations,  the  incident  wave  is  assumed  to  un¬ 
dergo  a  phase  shift  of  less  than  tt  in  traversing  the  scatterer. 
In  general,  however,  the  scatterer  cannot  be  assumed  to  be 
weak.  Also,  in  practice,  scattering  is  frequently  in  the  reso¬ 
nance  region  in  which  the  wavelength  is  less  than  or  compa¬ 
rable  to  the  characteristic  size  of  the  scatterer.1  Under  such 
conditions,  ||7]|oo  cannot  be  considered  to  be  small,  and  in¬ 
verse  problems  in  the  resonance  region  are  improperly  posed 
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and  nonlinear.  The  multiple  scattering  interactions  between 
the  propagating  wave  and  the  inhomogeneity  must  be  taken 
into  account.  Linearization  is,  therefore,  not  possible,  and  the 
inverse  solutions  are  usually  obtained  by  iteratively  minimiz¬ 
ing  a  suitably  constructed,  regularized  objective  function. 

The  scattered  fields  measured  at  the  detectors  are  always 
corrupted  by  noise,  the  sources  of  which  are  physical,  but 
can  also  arise  from  the  numerics.  The  reconstructed  image  is, 
therefore,  also  noisy.  Moreover,  the  propagation  of  noise  de¬ 
pends  upon  the  mathematical  procedure  by  which  the  data 
are  transformed  into  the  image,  that  is,  the  reconstruction 
algorithm.  Thus  the  understanding  of  how  noise  in  the  data 
propagates  to  the  image  via  a  given  reconstruction  algorithm 
is  clearly  an  important  problem,  not  only  for  understanding 
the  characteristics  of  the  noise  in  the  reconstructed  images 
produced  by  the  algorithm,  but  also  for  comparing  the  per¬ 
formances  of  reconstruction  algorithms  for  given  noise  char¬ 
acteristics  in  the  measurements.  The  propagation  of  noise 
under  the  weak  scattering  condition,  where  the  linearized 
Born  or  Rytov  approximation  holds,  has  been  investigated  by 
several  researchers. 8-11  In  this  paper,  the  attention  is  focused 
on  the  nonlinear  inverse  scattering  of  an  arbitrary  scatterer 
with  frequencies  in  the  resonance  region. 

The  paper  is  organized  as  follows.  The  direct  problem 
and  its  Frechet  differentiability  are  briefly  reviewed  in  Sec. 
II.  The  inversion  algorithm  and  the  calculation  of  the  gradi¬ 
ent  of  a  nonlinear  functional  via  the  adjoint  fieldl2~U  are 
presented  in  Sec.  III.  The  inverse  scattering  covariance  ma¬ 
trix  of  the  image  noise  is  derived  in  Sec.  IV,  and  its  relation 
to  that  of  the  linear  x-ray  computed  tomography  appears  in 
Sec.  V.  In  Sec.  VI,  the  inverse  scattering  covariance  matrix  is 
specialized  to  the  limiting  case  of  the  Bom  approximation. 
Section  VI  also  includes  a  short  discussion  about  the  appli¬ 
cability  of  the  present  analysis  to  obstacle  problems  and 
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physical  optics  approximation.  Finally,  a  brief  summary  of 
the  work  is  presented  in  Sec.  VII. 


II.  THE  DIRECT  PROBLEM  AND  FRECHET 
DIFFERENTIABILITY 

The  scattering  is  described  by  the  Lippmann-Schwinger 
integral  equation, 1-4  which  is 

^(x)  =  r(x)-k20f  G(0)(*|; y)  700000,  Cl  e  Rd, 

■In 

(1) 

d  =  2,3, 

subject  to  Sommerfeld’s  radiation  condition : 

lim|x|_,00|x|m[<?|,c|^c  -  ik0iffc]  =  0, 

ni  =  {d- 1)/2,  and  the  limit  is  uniform  in  |jc|  in  all  directions, 
where  d  is  the  dimensionality  of  the  problem.  17  is  the  com¬ 
pact  support  of  y,  that  is,  a  bounded  region  in  space.  fi"(x) 
=elk o%*  k0ESd~l,  the  unit  sphere  in  Rd,  is  the  incident 
plane  wave.  ij/=i//n+i/^c  is  the  total  field,  if/0  being  the  scat¬ 
tered  component.  From  Eq.  (1). 

il?c(x)  =  -kl  \  G(Q\x\y)yiy)ip{y),  x  $  17.  (2) 

■In 

In  Eqs.  (1)  and  (2),  y(x)  =  1  -(c0/c(x))2  describes  the  spatial 

variation  of  the  unknown  inhomogeneity  to  be  recovered, 

and  the  support  of  y  is  assumed  to  be  compact.  c(x)  is  the 

speed  of  sound  in  the  object.  G^(x|y)  is  the  continuous 

free-space  Green’s  function.  G^(x|y)  =  (;74)//[)1)(A:o|x-y|) 

in  two  dimensions  (2D),  and  G(0)(.r| y)={\!  ATr)elk^x~y^!i\x 

—  v I )  in  three  dimensions.  h\}]  is  the  zeroth-order  Hankel 

r  15  ^ 

function. 

The  solution  of  Eq.  (1)  can  be  written  as 

fi=£~xifi".  (3) 

In  Eq.  (3),  £  =  I+Qy.  L2(I7)— >L2((7)  is  the  Lippmann- 
Schwinger  operator.  I  is  the  identity,  and  Q  the  Green  func¬ 
tion  operator.  Furthermore,  Eq.  (2)  can  be  recast  as 

F(y)  =  i/fc.  (4) 

The  scattering  operator  F(y)  in  Eq.  (4)  is  Frechet  differ¬ 
entiable.  That  is,  there  exists  a  linear  operator,  F' ,  the 
Frechet  derivative  of  F,  such  that  \\F{y+h)-F(y)-F'h\\ 
~o(||/z||),  h  eL2(CI)  being  an  arbitrary  vector.  Some  details 
about  the  existence  and  uniqueness  of  the  Frechet  derivative 
can  be  found  in  Hohage.16 

The  mathematical  form  of  the  derivative  is  of  impor¬ 
tance  to  the  discussions  that  follow.  A  terse  derivation  is, 
therefore,  presented.  Let  fi=yijj.  Then  F{y)  =  Qijj,  and  Eq.  (1) 
becomes 

=  (5) 

in  which  £'  =/-  yQ  is  the  operator  adjoint  to  Also, 


tfin  +  gfi,  (6) 

Differentiating  F(y)  =  (jij/  formally  in  y  yields:  F’  =  Qi//' . 
From  Eqs.  (5)  and  (6),  it  is  not  difficult  to  show  that 
£\fi' h)  =  if/h,  from  which  it  follows  that  F'  =g[£l]~x  if/. 

In  finite  space  dimensions,  the  operators  Q,  C,  and  £ 
are  replaced  by  their  corresponding  matrix  operators,  which 
are  denoted  here  by  G,  L,  and  L+,  respectively.  We  then  ob¬ 
tain  the  Jacobian  matrix  of  the  scattered  field,  namely 

/sc  =  G[L+]-1A^,  (7) 

and  its  transpose 

=  (8) 

A  is  used  to  denote  a  diagonal  matrix.  Thus  {A^}i;  = 

Let  M  be  the  number  of  the  detectors,  and  the  grid  be  of  size 
NXN.  Then  Jsc  is  a  MXN2  matrix. 

III.  THE  INVERSION  ALGORITHM 

The  solution  of  the  inverse  scattering  problem  was  ob¬ 
tained  by  minimizing  the  following  nonquadratic,  Tikhonov 
type  functional: 

<F(r) = d/2  ^  -  re) + «iir-  y\\l 

(9) 

t[fne  is  the  measured  (noisy)  field,  2  the  covariance  matrix 
of  the  noisy  measurements,  and  a>0,  a  small,  the  regu¬ 
larization  parameter,  y  is  the  initial  estimate.  It  is  known17 
that  the  regularizing  term,  as  written  in  Eq.  (9),  contributes  to 
the  smoothness  of  the  solution.  It  is  to  be  noted  that  the 
functional,  <T>(  y),  is  nonquadratic  since  the  scattering  opera¬ 
tor,  F{y)  in  Eq.  (4)  is  nonlinear. 

Let  <J>  y=Vy<J>  be  the  functional  derivative  of  (b  with 
respect  to  y.  The  basic  reconstruction  is  then  given  by  the 
solution  which  is  the  null  point  of  the  gradient.  However,  in 
order  to  calculate  the  error  propagation  from  data  to  image,  it 
is  necessary  to  specify  the  mathematical  form  of  the  func¬ 
tional  derivative.  In  the  present  work,  the  derivative  was  ob- 

12 

tained  via  the  use  of  the  adjoint  fields.  The  result  is 

NW 

<^>,y=Re2  [A^.t/fJ]  +  2a(y—  y),  (10) 

i=  i  1 

where  Nw  is  the  total  number  of  the  waves  incident  on  the 
scatterer,  and  i/fj  is  the  adjoint  field,  if/l  is  given  by 

<A)-  =  L-Vjin,  i/ffn=G7j.  (11) 

In  Eq.  (11),  Tj  =  ij/f  - 1//)'6  is  the  complex  residual  vector  at 
the  detectors.  Its  mth  component,  -  <//™,  the  complex 
residual  at  the  mth  detector,  is 

rjm  =  [Re(r ;m)/Re(oym)]  +  ([Im(r7J/Im(cjjj] .  (12) 

Re(Im)  is  the  real  (imaginary)  part  of  a  complex  variable, 
and  the  overbar  implies  complex  conjugation. 

A  comparison  of  Eq.  (11)  with  the  Lippmann-Schwinger 
solution  in  Eq.  (3)  shows  that  t/f  can  be  obtained  from  Eq. 
(1)  if  if/jm  is  substituted  for  ifi"  in  that  equation.  Then  from 
Eq.  (10)  it  follows  that  the  computation  of  the  functional 
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derivative,  $  y,  of  the  objective  function,  $  in  Eq.  (9),  re¬ 
quires  the  solutions  of  only  two  forward  problems,  one  for  if/, 
and  another  for  if/1',  for  any  incidence.  It  should  be  mentioned 
at  this  point  that  a  closely  related  problem  occurs  in  nonlin¬ 
ear  tomography  where  the  derivative  of  a  nonlinear  operator 

is  calculated  via  the  solution  of  a  linear  adjoint  problem, 

18 

called  adjoint  dijferentiation. 

The  reconstruction  algorithm,  therefore,  consists  in  solv¬ 
ing  the  system:  <1>  r=0.  The  solution  of  the  inverse  problem 
then  reduces  to  that  of  locating  the  zeros  of  a  nonlinear  sys¬ 
tem  of  equations  given  in  Eq.  (10).  It  is  this  algorithm  that 
forms  the  basis  for  the  analysis  of  the  image  noise  in  this 
work. 


IV.  THE  PROPAGATION  OF  NOISE 

Let  £=(£)  + be  a  random  quantity  in  which  (ij)  is  the 
ensemble  mean  and  the  noise  component.  Following  this 
notation,  y={y)  +  ey,  if/=(if/)+ e*,  ^sc=(^sc)  +  esc,  ^  =  (1^) 
+  et,  and  i//ne=(tf/ne)+n.  is  the  noise-free  data,  and 

hence  (i/'me)  =  </£c,  if/fj  denoting  the  field  scattered  by  the  true 
object.  Implicit  is  the  assumption  that  the  magnitudes  of  the 
fluctuations  are  small  compared  to  their  mean  values.  Equa¬ 
tion  (10)  now  becomes 

$r=Re[<^)  +  <^et  +  ^t)eM] 

+  [higher  order  terms]  +  2a((y)  -  y  )  +  2 aey. 

(13) 

For  the  simplicity  of  notation,  the  sums  over  the  index  j  for 
the  incident  field  were  omitted  in  Eq.  (13).  These  will  be 
introduced  at  the  end.  Setting  the  right-hand  side  of  Eq.  (13) 
to  zero,  and  separating  the  mean  and  the  random  part,  yields 

Re[<'AX'At)]  +  2a((y)  -  y  )  =  0, 

for  the  ensemble  mean,  and 

+  {if/)e'lr\  +  2aey=Q  (14) 

for  the  random  part.  We  also  introduce  the  assumption  that 
the  iterations  have  converged  to  within  a  S  neighborhood  of 
the  solution,  in  which  case,  the  term  {if/)e^  in  Eq.  (14)  can 
be  neglected,  being  of  the  higher  order.  Equation  (14)  then 
reduces  to 

Re[Awe+]  +  2aAy=0,  (15) 

where  A,  as  earlier,  denotes  a  diagonal  matrix. 

Noting  that  ( i//1")  —  0  by  virtue  of  the  above-introduced 
assumption,  the  random  part,  e\  of  the  adjoint  field,  <//,  is 
given  by  the  relation:  et=(L)_l  efiri,  where  is  the  random 

part  of  (//"’= GF,  as  defined  in  Eq.  (11).  Moreover,  (L)  is  the 
Lippmann-Schwinger  operator  in  which  y=(y),  i.e.,  (L)  =  I 
-GA(Yy  The  quantity,  F,  defined  in  Eq.  (12),  contains  both 
the  data  noise,  n,  as  well  as  the  noisy  component  of  the 
scattered  field.  Let  F  be  written  as:  F=Fc+n,  where  Fc  de¬ 
notes  the  latter.  Then 


ef  =  (L)-1G[Fc  +  n],  (16) 

where  Fc  is  the  noise  component  of  the  scattered  part  of  F 
and  n  is  the  data  noise.  From  the  definition  in  Eq.  (12)  for 
F: 


= 

jm 


Re  e? 


jm 


Ime 


MIL 


Re  ojm  1m  ojm 


(17) 


and  a  similar  expression  applies  to  n.  From  Eqs.  (8)  and  (16), 
one  obtains 


Re[AwF]  =  R  etJ^  +  n)].  (18) 

In  order  to  proceed  further,  the  scattered  esc  is  expressed 
in  terms  of  the  object  ey.  It  is  readily  verified  that 
=Jscey.  In  view  of  Eq.  (2), 

F0  =  G(  y)  e^.  (19) 

Now  replacing  y=(y)  +  ey  and  if/=(if/)  +  e^  in  the  Lippmann- 
Schwinger  Eq.  (1),  separating  out  the  ensemble  averaged  and 
the  random  part  on  both  sides  of  the  result,  it  is  readily 
obtained  that 

6*=\{L)-lGA^\6y  (20) 

to  the  first  order.  Substituting  Eq.  (20)  in  Eq.  (19)  then  yields 
eiC=Jscey,  upon  using  Eq.  (7).  The  identity  also  follows 
from  the  definition  of  the  Jacobian  matrix. 

We  next  use  the  relation  esc=7scey,  in  Eq.  (17),  and  re¬ 
place  the  result  in  Eq.  (18).  Substituting  the  outcome  in  Eq. 
(15),  and  upon  setting  the  result  to  zero  yields 

(Re(/src)Re(2-2)Re(/sc)  +  Im(/[c)Im(£-2)Im(./sc) 

+  2  a I)  ey  =  Re(7fc)Re(X_2)Re(F) 

+  Im(./^c)Im(2_2)Im(«) .  (21) 

In  Eq.  (21),  2  denotes  the  covariance  matrix  of  the  data 
noise.  Thus  X;y-=oy,-. 

In  order  to  simplify  the  notations,  let  us  define  P,[J] 
^■{[S-'I/sc]},  and  1=1,2.  The 

subscripts  on  P  are  used  to  indicate  the  real  and  imaginary 
parts  of  the  operator.  Thus  P;=1(P,= 2)  imply  that  only  the  real 
(imaginary)  parts  of  the  operators  in  the  argument  of  P  are  to 
be  considered.  Moreover,  we  now  sum  over  the  number  of 
incident  waves  Nw  in  the  above-presented  expressions,  and 
define  A=Xf='yX2=1P,(^]P,{iJ],  and  P=X^X2=1PI(/7]PI(nJ.]. 
Equation  (21)  then  takes  the  form:  [A  +  2aI\ey=Bn,  from 
which 

ey=[A  +  2aI]~lBn.  (22) 

Now  cov(y)  =  (er{er}7).  From  Eq.  (22),  it  follows  immedi¬ 
ately  that 

cov(y)  =  [A  +  2aI]-1(Bn{Bn}I)[[A  +  2  a/]7]"1 .  (23) 

For  a  nonstochastic  scatterer  such  as  an  acoustic  inho¬ 
mogeneity,  it  is  reasonable  to  assume  that  data  are  uncorre¬ 
lated  with  respect  to  the  receiving  transducers  and  view 
angles.  The  scattering  data  are,  however,  complex,  and  their 
real  and  imaginary  parts  are  both  corrupted  by  noise.  This  is 
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(28) 


expressed  through  the  complex  cr.  A  correlation  can  exist 
between  crRe  and  <xIm.  However,  for  a  nonrandom  acoustical 
scatterer,  it  is  reasonable  to  neglect  this  correlation  (see  also 
Ref.  11).  In  other  words,  co v(njm,n,jimi)  =  o28jjtSmmi.  The  co- 
variance  matrix  in  Eq.  (23)  then  reduces  to 

cov(y)  =  [A  +  2a/]"1A[[A  +  2a/]-1]7'.  (24) 

Equation  (24)  is  the  covariance  matrix  of  the  image  noise  in 
the  inverse  scattering  algorithm  given  in  Eq.  (10). 


V.  LINEAR  PROBLEM  AND  REAL  DATA:  X-RAY  CT 

The  P  operators  in  Sec.  IV  are  nonlinear,  i.e.,  they  de¬ 
pend  nonlinearly  on  the  solution  (y).  This  is  explicitly  indi¬ 
cated  in  the  fundamental  equation  of  scattering,  namely, 
i l/sc=F(y),  [Eq.  (4),  Sec.  II].  F(y)  indicates  that  F  is  a  non¬ 
linear  functional  of  y.  Contrarily,  if  the  operator  F  acts  lin¬ 
early  on  y,  as  for  example,  if  F  is  a  matrix  acting  on  a  vector 
y,  then  it  is  written  as  Fy.  In  Eqs.  (7)  and  (8)  for  the  Jaco- 
bians,  the  nonlinearity  is  expressed  through  and  (L), 
both  of  which  are  (y)-dependent.  The  complex  nature  of  the 
scattered  data  introduces  a  further  layer  of  complication. 
Clearly,  the  covariance  matrix  must  simplify  if  the  operator 
is  linear  and  the  data  g  are  real,  that  is,  g  =  Fy  instead  of  g 
=F(y)  as  in  Eq.  (4). 

An  important  case  of  a  linear  problem  with  real  data  is 
that  of  the  x-ray  computed  tomography.19'"0  In  x-ray  CT,  the 
scattering  matrix  F  becomes  the  projection  matrix  //,  the  i/th 
element  of  which  represents  the  intersection  of  the  jth  x  ray 
with  the  ith  resolution  element.  The  noisy  x-ray  CT  data  can 
be  written  as  g=Hy+n=(g)+n.  In  this  case,  Eq.  (24)  for 
cov(y)  simplifies  to 

cov(y)  =  [HTt~2H+  2a]-1[//rS-2//] 

X  \[HT1r2H  +  2a]"1]7]  (25) 

A  special  case  of  Eq.  (25)  is  the  so-called  MAP  ( maxi¬ 
mum  a  posteriori )  x-ray  CT  reconstruction  algorithm21  in 
maximum  entropy  with  Poisson  noise  in  the  data.  The  corre¬ 
sponding  objective  function  is  given  by 

<b(y)  =  log{/?(g|y)}  +  AK(y),  (26) 


[//rX"2//  +  a2R\ey  =  HTl-2n . 

From  Eq.  (28),  the  covariance  for  x-ray  CT  follows  immedi¬ 
ately.  It  is  to  be  noted  that  for  Poisson  noise,  2  is  the  same  as 

A<s>- 

A  comparison  of  Eqs.  (24)  and  (28)  shows  that  the  error, 
ey,  has  the  same  mathematical  form  in  both  x-ray  CT  and 
acoustic  inverse  scattering,  namely, 

ey=[C  +  R(a)]-lEn, 

where  R(a)  is  the  regularizing  term.  However,  in  the  case  of 
linear  x-ray  CT,  C  and  E  involve  the  linear  operator  of  pro¬ 
jection,  H,  and  its  transpose.  In  acoustic  inverse  scattering, 
however,  C  and  E  are  relatively  more  involved,  being  non¬ 
linearly  dependent  in  y.  Further  details  and  an  application  of 
Eq.  (28)  to  a  dynamic  SPECT  problem  appear  in  Ref.  22, 
and  some  related  results  can  be  found  in  Refs.  23-25. 


VI.  THE  LIMIT  OF  WEAK  SCATTERING 


In  this  section,  Eq.  (24)  for  the  inverse  scattering  cova¬ 
riance  matrix  is  analyzed  in  the  limit  of  weak  scattering  in 
which  the  linearized  Born  approximation  holds.  Let  us  first 
assume  that  Re  cr/m=Im  a)m=  cr,  V/,  m,  and  therefore, 
Re  22=Im  2 2=cr2I.  22  is  a  diagonal  matrix:  (22)„  =  of.  From 
now  onward,  attention  will  be  confined  to  two  dimensions. 

Let 

Nw 

J=  2  .  (29) 

/=i 


In  the  Born  approximation,  Eq.  (8)  reduces  to 


J 


"A jG. 


(30) 


mA;  is  the  diagonal  matrix  for  the  jth  incidence.  That  is, 
{mAj)kk=exp{ik0dj-xk),  dj  being  the  unit  vector  specifying 
the  direction  of  the  jth  incident  wave.  Also,  the  far- field 
Green’s  function  has  the  form: 

-‘Wt  r  ia  _ 

— e~‘ lr/4  \  trk0Jl  ( k0a ) 


Jkt  - 


\R 


(31) 


in  which  use  was  made  of  the  asymptotic  expression:1 

//[0(f)  =  ^-ei(t-nirl2~vlA)\  i  +  (j)  |,  °°  . 


where  IZ(y)  is  a  positive  definite  symmetric  matrix  regular- 
izer.  p(g  |  y)  is  the  likelihood  of  data  g  given  a  distribution  y, 
assumed  to  be  positive.  For  Poisson  noise, 


P(g\y)  =  n^'|  [exp^Hj  Hj{y)j 


VhHuyd* 


gp 


(27) 


Note  that  the  quantity,  gp  in  the  denominator  in  Eq.  (27)  is 
an  integer  number,  being  the  number  of  counts.  MN  is  the 
total  number  of  data,  N  being  the  size  of  the  computational 
domain,  and  M  the  number  of  view  angles.  Carrying  out  the 
calculations  of  <J>  y,  with  the  log-likelihood  function , 
log  p{g |  y),  pig |  y)  in  Eq.  (27),  (for  details,  see  Ref.  22),  it 
is  obtained  that 


In  Eq.  (31),  d(  is  the  unit  vector  in  the  direction  of  the  €th 
detector  and  R  is  the  radius  of  the  detector  ring.  The  compu¬ 
tational  domain  was  assumed  to  have  been  discretized  into 
N2  circular  resolution  elements  following  the  widely  used 
discretization  scheme  of  Richmond."6  The  quantity,  a,  in  Eq. 
(31),  is  the  radius  of  such  a  circular  pixel  in  the  present,  2D, 
case.  The  geometry  of  the  so  discretized  computational  do¬ 
main  is  shown  in  Fig.  1.  In  Fig.  1,  xk  is  the  coordinate  of  the 
center  of  the  Zth  circle  in  the  grid,  and  dj,  df  are  the  unit 
vectors  along  the  jth  incident  wave  and  the  €th  detector, 
respectively.  The  detectors  are  located  on  a  ring  of  radius  R 
surrounding  the  scatterer. 

From  Eqs.  (30)  and  (31),  we  obtain  the  Jacobians  in  the 
Born  approximation,  namely 
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FIG.  1.  A  schematic  of  the  geometry  for  the  Bom  calculation  illustrating  the 
discretization  of  the  computation  grid  into  circular  region.  xk  is  the  location 
of  the  center  of  the  &th  resolution  element  in  the  grid,  dp  d(  are  the  unit 
vectors  along  the  jth  incident  wave  and  the  €th  detector,  respectively.  The 
detectors  are  located  on  a  ring  of  radius  R  surrounding  the  scatterer.  In  the 
numerical  computation  in  this  work,  the  grid  was  discretized  into  32  X  32 
circular  elements  each  of  radius  a.  R/ a  was  200.  That  is,  the  detectors  were 
placed  on  a  radius  which  was  200  times  the  radius  of  a  circular  pixel. 


eW 

[jlj]km=^R 


Y e  ,m'[TTk0Ji(koa) 


,ik0(drdjxk 


(32) 


and 


[Jj]mt=lR. 


■  '^emn'jTTk0Jl(k0a) 


e-ik0(drdm)xe 


Replacing  Eqs.  (32)  and  (33)  in  Eq.  (29),  and  after  some 
straightforward  algebra,  it  follows  that 


Jkt-C 


Nw  M 

ReX  X  e~ik°^di~dm,  (-Xk~Xt) 

j=  i  m=l 


(34) 


where  the  constant  C  =  (Trk0a2 / 2R)[J ^IcqCi)]2 .  Note  that  Eq. 
(34)  reflects  the  well-known  fact  that  in  the  Born  approxi¬ 
mation  the  Fourier  frequencies  of  the  object  are  confined 
within  a  circle  (in  2D)  of  radius  2 k0,  and  therefore,  the  maxi¬ 
mum  frequency  that  can  be  present  in  the  recovered  object  is 
bounded  above  by  2k0. 

Typically  in  practice,  the  product  NWM  is  a  large  num¬ 
ber  for  any  realistic  discretization  of  the  computational  do¬ 
main.  For  example,  for  a  64  X  64  size  of  the  computational 
grid,  this  product  is  4096  at  the  minimum.  In  practical  com¬ 
putations,  this  number  is  often  larger  as  the  system  is  almost 
always  overdetermined  in  actual  numerics.  We,  therefore,  ap¬ 
proximate  the  sums  in  Eq.  (34)  by  integrals.  Moreover,  the 
sums  over  j  and  m  are  independent  sums.  Assuming  without 
any  loss  in  generality  that  NW=M,  each  individual  sum  can 
be  approximated  by 


FIG.  2.  (a)  The  1024  X  1024  S  matrix  [Eq.  (36)]  for  a  32  X  32  grid,  (b)  The 
corresponding  T  matrix  with  a=  10-4.  The  block  Toeplitz  block  structure  of 
the  matrices  are  visible. 


Nw 

X  cos(z  COS  (j) )  ~ 
7=1 


cos(z  cos  4>)d<t>  =  NwJ0(z) , 


(35) 


recalling15  that  the  integral  in  Eq.  (35)  is  one  of  the  repre¬ 
sentations  of  the  cylindrical  Bessel  function  J0.  Also,  in  Eq. 
(35),  z=kQ\xk-x(\. 

In  view  of  Eq.  (35),  Eq.  (24)  for  the  covariance  reduces 
to  the  following  expression  in  the  Born  approximation, 
namely 

[covB(y)j  =  cr[7~1S7~1] .  (36) 


In  Eq.  (36),  S=CN2WZ,  ZH=  jl(k0\xk-Xf  |),  and  the  constant 
C  was  defined  earlier.  Moreover,  T=S+2aI.  Note  that  S  and 
T  are  each  a  block  Toeplitz  Toeplitz  block  or  a  BTTB 
matrix.27 

It  should  be  mentioned  that  the  object  enters  the  calcu¬ 
lation  through  the  noise  in  the  data.  In  x-ray  CT,  for  example, 
the  object  information  is  contained  in  S-2  in  Eq.  (28).  Simi¬ 
larly,  in  the  Born  approximation,  under  the  assumptions 
made  in  the  derivation  of  Eq.  (24)  (essentially,  the  same  as¬ 
sumptions  also  appear  in  Ref.  1 1),  the  object  appears  through 
cr  in  Eq.  (36)  which  may  originate,  say,  from  speckles.  Thus 
the  behavior  of  the  covariance  is  determined  completely  by 
Eq.  (36).  It  is  illustrated  here  assuming  a  32X32  computa¬ 
tional  grid. 

The  BTTB  S  matrix  (1024X  1024)  is  shown  in  Fig.  2(a), 
and  the  corresponding  BTTB  T  matrix  (1024X1024)  is 
shown  in  Fig.  2(b)  with  the  regularizing  parameter  a=10-4. 
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FIG.  3.  The  behavior  of  the  covariance  matrix,  cov(y)  [Eq.  (36)]  with  pixel 
separation.  The  variance  or  the  diagonal  elements  of  the  matrix  are  constant. 
Panel  (a)  corresponds  to  the  case  in  which  \/a= 2,  whereas  (b)  is  for  \/a 
=  3.  Covariance  is  peaked  at  the  center  and  falls  off  rapidly  with  the  inter¬ 
pixel  separation. 


Toeplitz  structures  are  clearly  visible  in  Fig.  2.  The  covari¬ 
ance,  covB(y),  given  in  Eq.  (36),  is  shown  in  Fig.  2.  The 
variance  is  found  to  be  constant  in  each  pixel,  whereas  the 
off-diagonal  elements  of  the  covariance  matrix  are  highly 
peaked  around  the  center  pixel,  and  decay  rapidly  in  an  os¬ 
cillatory  manner  as  the  interpixel  separation  increases.  In 
Fig.  3(a),  X/a  was  2,  while  in  Fig.  3(b),  A.  was  three  times  the 
radius  of  a  resolution  element.  In  a  recently  published 
report,11  the  noise  properties  of  Born  objects  were  calculated 
using  direct  Born  inversions.  It  is  known'  that  Born  data 
can  be  inferred  by  filtered  backpropagation  (FBPP),  which  is 
the  analog  of  filtered  backprojection  of  x-ray  CT19"0  in 
acoustical  scattering.  The  behavior  of  the  Born  covariance 
matrix  in  Ref.  1 1  using  FBPP  inversion  is  qualitatively  simi¬ 
lar  to  that  in  Fig.  3.  In  both  cases,  the  constant  diagonal  (the 
variance)  and  rapidly  decaying  oscillatory  off-diagonal  ele¬ 
ments  were  obtained.  Moreover,  both  methods  produced  cr 
dependence  in  the  covariance  matrix.  However,  the  depen¬ 
dence  on  k0  is  relatively  more  indirect  in  inverse  scattering 
than  in  FBPP  which  produced  a  simple  dependence.  Since 
the  inversion  algorithms  are  different,  only  a  qualitative 
agreement  is  to  be  expected. 

Remark.  The  preceding  analysis  of  the  propagation  of 
noise  from  the  data  to  the  reconstructed  image  was  carried 
out  in  the  framework  of  the  Lippmann-Schwinger  integral 
equation  of  scattering.  The  Lippmann-Schwinger  equation 
provides  the  appropriate  framework  for  scattering  calcula¬ 
tions  from  inhomogeneities,  that  is,  when  the  scattering  po¬ 
tential  has  no  discontinuity  across  its  support.  However,  the 
analysis  can  also  be  extended  to  the  obstacle  problem  in 
which  the  potential  is  discontinuous  across  its  boundary.  For 
the  obstacle  problem,  the  appropriate  framework  of  analysis 


1  3 

is  the  Helmholtz  representation,  ’  which  involves  surface  in¬ 
tegrals.  However,  the  inversion  can  be  formulated  into  a 
problem  in  nonlinear  parameter  estimation1  involving  objec¬ 
tive  functions  which  are  essentially  similar  to  $  in  Eq.  (9). 
The  parameters  to  be  estimated  in  the  case  of  the  obstacles 
are  the  variables  which  parametrize  the  surface,  e.g.,  spheri¬ 
cal  harmonics.  The  algorithm,  Eq.  (10),  remains  unchanged, 
and  the  functional  derivative  (now  in  the  parametric  vari¬ 
ables  of  the  surface)  of  the  objective  functional  can  again  be 
obtained,  as  in  the  case  of  the  inhomogeneities,  via  the 
method  of  adjoint  fields,  and  again  requiring  the  solutions  of 
only  two  forward  problems.  The  details  of  the  procedure 
appear  in  Ref.  28.  However,  in  the  case  of  obstacle  scatter¬ 
ing,  the  derivatives  are  not  the  regular  partial  derivatives  (as 
in  Sec.  Ill),  but  are  the  so-called  shape  derivatives ,29~33 
which  are  defined  in  the  following  manner.  If  V:R3— >R3  is  a 
vector  field  that  deforms  a  domain  O0  e  /?3— e  R  \  t  be¬ 
ing  the  perturbation,  then  the  first-order  shape  derivative  of 
4>,  #1],  is  defined  by 


i =  limit, 


t{i  +  tv){x)-ux) 

t 


VXe  R 3. 


if/0,  if/,  denote  the  fields  for  the  unperturbed  and  perturbed 
boundary,  respectively,  corresponding  to  the  boundary  con¬ 
dition  of  the  scattering  problem,  if)  1  ^  is  the  total,  material,  or 
substantial  derivative.  Moreover,  </<W  =  ip'  +  V-  V<p,  i f/'  being 
the  partial  derivative,  namely,  ijt'  =limit,^0[iArW- <Ao(X)]/F 
In  a  similar  vein,  the  analysis  can  also  be  extended  to  the 
limit  of  the  physical  optics  approximation.  It  is  to  be  noted 
that  the  physical  optics  approximation  bears  the  same  rela¬ 
tion  to  the  Helmholtz  representation  as  the  Born  approxima¬ 
tion  does  to  the  Lippmann-Schwinger  integral  equation.  In 
both  cases,  the  total  field  within  the  surface  integrals  is  re¬ 
placed  by  the  incident  waves  weighted  by  suitable  constants. 


VII.  CONCLUSIONS 

The  propagation  of  noise  from  the  data  to  the  speed  of 
sound  image,  reconstructed  by  inverse  scattering  within  the 
framework  of  the  Lippmann-Schwinger  integral  equation  of 
scattering,  was  investigated.  The  inversion  algorithm  consti¬ 
tuted  in  locating  the  zeros  of  the  gradient  of  a  Tikhonov-type 
functional  in  the  unknown  speed  of  sound.  The  gradient  of 
the  functional  was  computed  by  the  method  of  the  adjoint 
fields.  An  analytical  expression  was  obtained  for  the  image 
covariance  matrix.  When  applied  to  the  case  of  a  linear  map¬ 
ping  with  real  data,  the  inverse  scattering  expression  reduced 
to  that  of  the  linear  x-ray  computed  tomography.  This  was 
demonstrated  by  considering  the  maximum  a  posteriori  like¬ 
lihood  algorithm  using  the  maximum  entropy  functional.  The 
full  inverse  scattering  covariance  was  also  investigated  in  the 
limiting  case  of  the  linearized  Born  approximation.  It  was 
found  that  in  this  linearized  approximation,  the  variance  re¬ 
mained  uniform  over  the  image,  whereas  the  off-diagonal 
elements  of  the  covariance  matrix  exhibited  a  fast  decaying 
oscillatory  behavior.  Similar  behaviors  were  also  obtained  in 
a  recent  report  on  Born  inversion  using  the  filtered  back- 
projection  algorithm.  Finally,  the  applicability  of  the  analysis 
presented  to  the  problem  of  obstacle  scattering  and  in  the 
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limit  of  the  physical  optics  approximation  was  briefly  dis¬ 
cussed. 
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Abstract 

The  propagation  of  acoustic  waves  in  a  circularly  cylindrical  waveguide  in  the  pres¬ 
ence  of  flow  was  investigated.  An  integral  equation  for  the  total  pressure  in  the 
duct  was  derived,  and  subsequently  reduced  in  the  Born  limit,  i.e. ,  for  flow  speeds 
much  slower  than  the  speed  of  sound  in  the  quiescent  medium.  The  Born  term  was 
then  used  to  determine  the  mode  phase  speeds,  and  hence  the  effects  of  flow  on 
wave  propagation  in  the  duct.  Results  are  reported  for  the  parabolic  laminar  and  a 
mixed  parabolic-flat  turbulent  velocity  distribution  of  a  homogeneous  fluid  flowing 
inside  ducts  of  Dirichlet,  Neumann  and  Robin  surfaces.  For  the  Neumann  bound¬ 
ary  and  mixed  flow,  the  integral  equation  yielded  results  similar  to  those  obtained 
by  recently  reported  perturbative  calculations.  However,  the  mathematics  is  more 
transparent  in  the  integral  equation  method  which  is  also  more  versatile,  and  can 
be  applied  with  relative  ease  to  numerous  practical  situations  which  may  involve 
inhomogeneous  media,  arbitrary  flow  profiles  and  boundaries. 
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1  Introduction 


The  propagation  of  acoustic  waves  in  circular  ducts  carrying  flow  is  an  impor¬ 
tant  problem  with  a  well-developed  literature  [1-13].  Various  techniques  exist 
for  the  solution  of  the  problem  including  series  solutions  of  the  governing  dif¬ 
ferential  equations  (e.g.,  [8,10]),  finite  element  calculations  (e.g.,  [6,12]),  and 
hydrodynamical  approach  (see  e.g.,  [9,13]).  Recently,  a  perturbative  approach 
was  reported  [11]  for  calculating  the  changes  in  the  mode  phase  speeds  due 
to  flow  in  a  sound-hard  duct  with  a  view  to  obtaining  accurate  flow  metering 
(see  also  [5]  for  this  application).  In  this  paper,  we  discuss  the  application  of 
the  integral  equation  approach,  well-known  [14, 15]  from  acoustical  scattering 
theory,  to  the  solution  of  the  wave  propagation  problem  in  a  duct  in  the  pres¬ 
ence  of  flow.  A  formal  integral  equation  is  derived  for  the  total  pressure  field 
in  the  flowing  fluid,  which  is  then  reduced  in  the  Born  limit.  In  this  limit,  it  is 
assumed  that  the  flow  speeds  are  much  slower  than  the  speed  of  sound  in  the 
quiescent  medium  (ie.  in  the  absence  of  flow).  The  method  was  used  in  order 
to  determine  the  mode  phase  speeds  in  the  presence  of  the  flow,  and  hence 
investigate  the  effects  of  the  flow  on  the  wave  propagation.  Results  were  ob¬ 
tained  for  the  laminar  parabolic  as  well  as  a  mixed  laminar-turbulent  velocity 
profile  of  the  fluid  flowing  in  ducts  with  a  Dirichlet  (sound-soft),  Neumann 
(sound-hard)  and  Robin  (impedance)  boundaries. 

In  the  Born  approximation,  with  Neumann  boundary,  the  integral  equation 
produced  results  similar  to  those  that  were  obtained  in  [11]  using  a  pertur¬ 
bative  approach.  The  mathematics  of  the  present  method  is,  however,  more 
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transparent.  Unlike  most  methods,  the  integral  equation  does  not  (a  priori ) 
assume  any  particular  form  of  the  pressure  held  except  its  time  harmonic- 
ity.  The  held  is  rather  a  result  of  the  calculations.  In  addition,  the  integral 
equation  formalism  is  applicable  to  problems  including  nonuniform  how  and 
arbitrary  velocity  prohles  and  boundary  conditions.  In  addition,  the  powerful 
techniques  that  exist  for  solving  various  acoustical  scattering  and  inverse  scat¬ 
tering  problems,  can  be  brought  to  bear  on  the  problem  of  how  under  quite 
general  conditions. 

The  plan  of  the  presentation  is  as  follows.  In  Section  2,  the  basic  equations  are 
presented,  and  a  formal  solution  for  the  total  pressure  held  in  a  flowing  huid  is 
developed.  The  solution  is  then  specialized  to  the  Born  limit.  Section  3  derives 
Green’s  functions  which  satisfy  the  boundary  condition  and  Sommerfeld’s  ra¬ 
diation  condition  at  infinity.  In  Section  4,  the  incident  pressure  and  the  Born 
approximated  total  pressure  held  are  derived.  Section  4  further  demonstrates 
how  the  Born  approximated  inhomogeneous  part  of  the  total  pressure  leads 
to  the  determinations  of  the  phase  speeds  of  the  propagating  modes  in  the 
presence  of  how.  Section  5  contains  the  details  of  the  numerics.  The  results 
obtained  are  then  discussed  in  the  following  Section  6.  The  text  is  concluded 
by  a  few  remarks  on  the  extension  of  the  integral  equation  approach  to  general 
how  prohles,  huid  media  and  waveguide  surfaces. 


2  The  Equations  of  Wave  Propagation 

Consider  the  how  of  a  homogeneous  huid  in  a  long  cylinder  of  radius  R.  Let 
V  =  {uj}f=1  describe  the  velocity  held  of  the  huid.  Then  t  =  t' ,  ay  =  ay  +  v{t' 
represent  the  coordinate  transformations  between  the  laboratory  (unprimed) 
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system  and  the  one  co- moving  with  the  flow  (primed).  In  the  primed  coordi¬ 
nate  system,  the  flow  equations  are: 


Up  »  U'L’i 

— — b  Po  /  ttt  —  0  :  momentum  conservation. 
at'  ^  ox\ 

dvi  dp 

Pq— - b  ttt  —  0  :  mass  conservation, 

at'  ox\ 


(1) 

(2) 


p  and  po  are  the  perturbed  and  unperturbed  mass  density,  respectively,  and 
p  is  the  pressure.  It  is  well-known  [16]  that  under  the  above  coordinate  trans¬ 
formation,  Eqs.  (1)  and  (2)  combine  to  yield: 


V2p  +  kl 


1  +  -V  ■  V 

u 


P  =  0, 


(3) 


assuming  time  harmonicity.  u  is  the  angular  frequency,  and  ko  =  oj/cq  the 
wavenumber,  Co  being  the  sound  speed  in  the  quiescent  medium.  In  the  first 
approximation  with  respect  to  the  Mach  number,  M  =  |Vmax|/c0,  and  as¬ 
suming  that  the  flow  is  quasi-steady  (typical  time  variation  r  >>  cu-1),  the 
second-order  term  in  Eq.  (3)  can  be  neglected  [17]  resulting  in  the  approxi¬ 
mation: 


V2p  +  kl 


2  i 

1  +  -V • V 

u 


p  =  0. 


(4) 


In  cylindrical  coordinates  r  =  {p,  0,  z },  V2  =  Vy  +  p~2d^  +  dzz  with  V^_  = 
dpp  +  {1/ p)dp.  It  is  further  assumed  that  the  flow  is  along  the  ^-direction 
which  is  also  the  directrix  of  the  cylindrical  waveguide  and  is  radially  varying. 
Therefore,  V  =  V(p,  0,  z)  =  zV  (p),  z  being  the  unit  vector  in  the  ^-direction. 
Equation  (4)  then  becomes: 


V2p  +  kl 


1  +  -V(p) 

IX 


^  =  o. 

oz 


(5) 
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Let  G^°)(r|r5)  be  the  “free-space”  Green’s  function,  that  is,  Green’s  function 
in  the  quiescent  fluid,  between  the  source  location,  rs,  and  the  observation 
point  r.  It  is  described  by: 

V2  +  kl]  G(°)(r|rs)  =  ~-8{p  -  ps)<5(0  -  </>s)6(; z  -  zs).  (6) 

Furthermore,  let  Gr^(r|r,s)  also  satisfy  the  boundary  condition  on  the  waveg¬ 
uide,  namely 

BGW  =  0  on  the  surface  of  the  cylinder,  (7) 

B  being  the  boundary  operator.  B  —  I,  I  the  identity  operator,  corresponds  to 
the  Dirichlet,  B  =  n- V,  n  the  unit  normal  on  the  boundary,  to  the  Neumann, 
and  B  =  n  ■  V  +  ia  to  the  Robin  boundary.  The  parameter  a  denotes  the 
quantity  R/3,  (3  =  (3'  +  i/3”  being  the  acoustic  admittance  of  the  cylinder  wall. 
In  terms  of  G^,  the  total  pressure,  p,  at  any  point  in  the  fluid,  can  be  obtained 
from  Eq.  (5)  as: 

p(r)=pmc(r)  +  /dr'  V(r^%lr')flp(r').  (8) 

pinc(r)  is  an  appropriate  incident  pressure  held,  and  is  the  volume  of  the 
waveguide.  Moreover,  in  Eq.  (8),  V(x)  =  —(2 ik0/c0)V(x). 

Equation  (8)  is  essentially  the  Lippmann- Schwinger  integral  equation  of  scat¬ 
tering  [14,15]  except  that  the  propagator  is  now  G^0^(r|r'){(9/9^}  instead  of 
G alone,  as  in  the  usual  scattering  equation.  It  is  noteworthy  that  the  how  it¬ 
self  (via  V(r))  provides  the  scattering  centers  and  acts  as  an  inhomogeneity,  no 
physical  particulate  matters  being  necessary  to  provide  the  centers,  although 
the  presence  of  additional  physical  inhomogeneities  (of  arbitrary  nature)  poses 
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no  problem  in  the  integral  equation  approach.  In  other  words,  the  flow  gives 
rise  to  what  may  be  termed  acoustic  velocity  inhomogeneity  for  diffraction  (see 
also  [18]).  Since  the  diffraction  scales  as  the  Mach  number  [19],  the  scattering 
from  the  velocity  inhomogeneity  will  indeed  be  small  for  M  <<  1,  and  the 
Born  approximation  can  be  invoked.  It  is  well-known  [14-16,20,21]  that  in  the 
Born  approximation,  Eq.  (8)  reduces  to: 

p( r)  =  pinc(r)  +  J  dr'  |f/(r')G'(0)(r|r')^:|  pinc( r').  (9) 

In  other  words,  the  Born  approximation  consists  in  replacing  the  total  pres¬ 
sure,  p,  in  the  integral  by  the  incident  field,  pinc.  The  integral  on  the  right-hand- 
side  of  Eq.  (9)  is  the  inhomogeneous  term  in  the  solution  for  the  total  pressure, 
p,  and  is  denoted  by  Pih  in  the  sequel.  Therefore,  p( r)  =  pinc(r)  +  Pih(r).  It  is 
Eq.  (9)  which  is  of  interest  here. 

The  free-space  Green’s  function,  G^°\  is  derived  next. 


3  The  ” Free-Space”  Green’s  Function 

As  was  already  pointed  out,  unlike  the  standard  Green’s  functions  in  an  un¬ 
bounded  space  ((i/4)HQ1\k0\x —  y |))  in  2-D,  H({)  1 }  being  the  cylindrical  Hankel 
function  of  order  zero  of  the  first  kind  [22]  and  exp(iko\x  —  y\)/4ir\x  —  y |)  in 
3D),  G (-0)  here  is  a  boundary  Green’s  function  which  satisfies  both  the  condi¬ 
tion  on  p  imposed  on  the  boundary  of  the  waveguide,  as  well  as  the  radiation 
condition  at  infinity.  For  a  boundary  of  arbitrary  shape,  this  Green’s  function 
may  be  difficult  to  obtain,  but  fortunately  for  boundaries  of  canonical  shapes 
(as  in  the  present  case),  it  can  be  obtained  rather  easily. 
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In  order  to  obtain  an  expression  for  G^°\  we  first  note  that: 


1  OO 

-  &)  =  5-  E 

■£=— OO 


(10) 


in  terms  of  the  definition  of  the  angular  Fourier  transform  pair,  (/(</>),  /(m))  : 

/«■)  =  E  /(my*,  /(m)  =  2-  f  f(4>)e-‘m*d4>. 

m=—oo  27F  4 

Moreover,  upon  using  the  orthogonality  relation  for  the  Bessel  functions: 

R 


J  d rrJn{Kisr)Jn(n'Csr)  =  0,  nts  ^ 


(11) 


which  is  known  to  hold  under  quite  general  conditions  [23],  the  radial  delta 
function,  (1/ p)5{p  —  ps ),  can  be  expressed  as: 


-&{p  ~  ps)  =  -Ly,  Mp)Mps), 

p  K  I. 


(12) 


L  in  Eq.  (12)  denotes  the  joint  indices,  {is}.  In  other  words,  J2l  =  )C<£o 
The  orthonormalized  Bessel  function,  Jl(z),  is: 


JL(z)  =  Al  Je  (  kl- )  , 


(13) 


the  normalization  constant,  Al,  being: 


„  R 


R 


dpp 


JfZ  |k4 


(14) 


The  eigenvalues,  kl,  in  Eqs.  (13)  and  (14),  are  obtained  from  the  operator 
equation, 

>C(kl)  =  0,  (15) 
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the  operator,  /C,  being  boundary  dependent  via  Eq.  (7).  From  Eqs.  (10),  (12), 
and  (6),  it  follows  that: 


v2  +  fcg]  G<°>(r|rs)  =  -  Jp  X  eu^-^gL(p,  ps)S(z  -  zs;  L),  (16) 

where  gL(p,  ps)  =  Jl(p)Jl(ps )■ 

Let 

9t(p,  Ps ;  z,  zs)  =  gL{p,  Ps)h(z,  zs),  (17) 

be  the  reduced  Green’s  function.  Equation  (16)  then  becomes: 

G(0)(r|rs)  =  -^J2eie(^s)9L(p,  Ps)9l(z,  zs).  (18) 

From  Eqs.  (18)  and  (16),  it  is  readily  seen  that  ()l{z,zs )  is  given  by  the 
equation: 

{dzz  +  A 2l}9l(z,zs)  = -5(z- zs),  (19) 

where  AL  =  \jk%~  (yy)  •  The  solution  of  Eq.  (19)  is: 

Si(2’Zs)  =  !x7<AiI*"'sI'  (20) 

Collecting  the  results  in  Eqs.  (17),  (18),  and  (20),  the  free-spaee  Green’s  func¬ 
tion  in  the  waveguide  is  finally  derived  to  be: 

G'»)(r|rs)  =  (21) 

■nR-  L  AL 


The  phase  speed  of  the  waves,  cl  =  cu/A l,  is  boundary  dependent  through 


4  The  Total  Pressure 


The  total  pressure,  p,  is  given  by  Eq.  (9).  The  incident  pressure,  pmc,  is  given 
by: 

pinc(r)  =Y.PL)eU*JL(p)eiKLZ  (22) 


corresponding  to  V  —  0.  Replacing  Eqs.  (21)  and  (23)  in  Eq.  (9),  and  after 
some  straightforward  algebra,  the  total  pressure,  p,  is  obtained: 


p{  r)  = 

L 


i  (q)Al' 


L  s' 


■  f  d z>ei^\‘-^\eiK’‘,% 


(23) 


In  Eq.  (23),  L'  =  {£,  s'},  and  L,  s'  =  {£,  s,  s'}.  Moreover,  the  quantity,  Cl,s'(R) 
is  given  by: 

Cl,AR)  =  —  /  d p'p'JL(p')V(p')JLI{p'). 
u  J 
0 


Let  us  next  define  the  quantity: 
2  i 


tl„. = -%pTR-h  H)/d2'e'All”vw' 

R 1 


(24) 


In  view  of  Eq.  (24),  the  inhomogeneous  part  of  the  solution,  p\\,  in  Eq.  (23) 
can  be  written  as: 


Pih(r)  =  EE  Rls'C’Ls'F Ls'  ■  (25) 

L  s' 

At  this  point,  let  us  note  that: 

G'0)*P™  =  Y.Y.Cl>p..‘,  (26) 

L  s' 
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in  which 

ci”i  =  jdppj'  (|>)  j,  (fP) , 

o 

and  *  denotes  the  operation  of  convolution.  Using  Eq.  (26)  in  Eq.  (25)  results 
in: 


Pih(r)  = 

L  s' 


cLyS, 

/ — »(o) 

L'L,s' 


(27) 


Since  from  the  orthogonality  relation  (11),  Cf^s,  =  (R2 /2)A2Ls,5ssi,  Eq.  (27) 
then  becomes: 


Pih(r, 


=  EA, 

L 


CL 

Cf 


(G(°) 


*  p 


J  L 


(28) 


Next  consider  a  single  term  in  the  summation  in  Eq.  (28).  In  other  words, 
consider  only  a  single  mode,  say,  L,  propagating  in  the  waveguide.  In  this 
particular  case,  we  have 

P1h(r)  =  A1l|y(G<»>»p1”)i.  (29) 


It  must  be  noted  that  Eq.  (29)  is  the  Born  approximation  of  the  unperturbed 
inhomogeneous  term,  A l(Cl/C^)p,  in  the  unperturbed  differential  equation 
of  propagation.  Therefore,  the  differential  equation  of  which  Eq.  (29)  is  the 
Born  approximation,  is: 


V2  +  kl  -  Dl  pL( r)  =  0, 


(30) 


in  which 


Dr  =  A, 


Cl_ 

cf 


(31) 
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From  Eq.  (30),  it  follows  that  the  phase  speed  for  the  mode  under  consideration 


in  the  presence  of  the  flow  is  cl  —  oo  /  y  —  DL.  The  contribution  to  the  phase 

speed  dne  to  the  flow  is,  therefore,  the  difference: 


(Ac)i  = 


(32) 


Equation  (32)  is  similar  to  the  result  which  was  obtained  in  [10]  using  a  per¬ 
turbative  approach. 


5  Numerical  Computations 

The  details  of  the  numerics  are  presented  in  this  Section.  Eigenvalues  were 
calculated  from  the  operator  Eq.  (15)  which  gives  Ji(kl)  =  0  for  the  Dirichlet, 
J[{kl)  =  0  for  the  Neumann  and  klJ[(kl)  ~  iaRJi(hiL)  =  0  for  the  Robin 
boundary.  Eigenvalues  for  the  Dirichlet  and  Neumann  cases  were  generated 
using  a  web  interface  [24]  to  Wolfram  Research  Inc.’s  Mathematica  software. 
For  the  Robin  boundary,  the  eigenvalues  were  computed  with  The  Mathwork’s 
MATLAB  software  using  Newton’s  method. 

The  corresponding  normalization  constants  are:  AL  =  J^+i(kl),  Al  = 
and  Al  =  J^kl)^  —  respectively. 

Two  values  of  the  quantity,  (  =  \/3\koR  namely,  0.3  and  1.2  are  reported  for 
the  impedance  boundary.  For  (  =  0.3,  (3  <  X/2ttR ,  and  for  (3  in  this  range, 
the  roots  of  the  Robin  boundary  can  be  obtained  analytically  (e.g.,  see  [16]): 


—  \Z‘2iRoi\  /cqo  —  fco  1  T-  'i’(3r)/ R,  (33) 


M 


N 


Kis  ~  K(.s 


1  - 


(/?"  +  i(3')k0R 
{<}2~s2  . 


.  A  R  _  a(0)  , 
>  Ives  —  ives  ' 


(F'  +  i^ko 


RA 


(o) 

is 


1  -  4r 


(34) 


Moreover, 


A(0)  - 

lyts  — 


\ 


/.'(-IV1 


kN' 

R 


(35) 


and  the  superscripts  R  (N)  indicates  Robin  (Neumann)  eigenvalues,  respec¬ 
tively.  The  impedance  eigenvalues  were  calculated  by  numerical  root  finding. 
For  C  =  0.3,  they  were  identical  to  those  computed  from  Eqs.  (33)  and  (34) 
thereby  verifying  the  accuracy  of  the  root  founding  algorithm.  However,  for 
C  =  1.2,  (3  was  not  small  compared  to  \/2irR ,  the  analytical  results  were  not 
valid,  and  the  numerical  roots  had  to  be  used. 

Following  [11],  two  flow  velocity  distributions  were  considered,  which  were: 


Vi(p)  =  2v  ^1  —  jpj  :  laminar  flow, 
vm(p)  =  avi(p)  +  (1  —  a)vt(p ) :  a  mixed  laminar  —  turbulent  flow. 


(36) 

(37) 


v  is  the  mean  flow  speed,  and  the  subscript  /,  m  denote  parabolic  laminar  and 
mixed  flow,  respectively.  Moreover,  again  as  in  [11],  a  was: 


a(Re ) 


(38) 


where  n  =  4,  and  Re  is  Reynold’s  number:  Re  =  2 pvR/p,  p,  p  being  the 
fluid’s  mass  density  and  viscosity,  respectively.  R=  0.006  m  in  the  numerics. 
Numerical  computations  were  performed  at  the  frequency  of  1  MHz. 
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6  Results  and  Discussions 


The  relative  changes,  (A c)l,  in  the  mode  phase  speeds  due  to  flow,  given  in 
Eq.  (32),  constituted  the  primary  results  in  this  work.  The  profiles  of  (Ac)^ 
vs.  v  are  shown  in  Figs.  1  through  6,  the  relevant  parameters  in  the  calcu¬ 
lations  described  in  the  figures.  Considering  the  large  amount  of  data  (two 
flow  profiles,  three  boundaries  and  a  large  kl  vector  with  £  =  0,  1  and  2,  and 
s  —  0, 1,2,  3, 4,  for  each  £)  involved,  a  limited,  yet  a  representative  set  of  re¬ 
sults  is  presented.  The  profiles  for  the  Dirichlet  boundary  are  plotted  in  Figs. 
1(a)  -  1(d),  and  the  Neumann  results  are  displayed  in  Figs.  2(a)  -  2(d).  The 
data  for  the  impedance  boundary  with  Q  =  0.3,  are  shown  in  Fig.  3,  whereas 
Fig.  4  shows  the  same  results,  but  for  (  =  1.2.  (3  is  real  in  Figs.  3  and  4, 
imaginary  (3  showing  qualitatively  similar  variations.  Moreover,  Figs.  3  and  4 
show  the  real  part  of  (A c)l- 


One  distinguishing  feature  of  the  plots  is  that  (Ac)l  (Re  ((Ac)l)  for  the 
impedance  condition)  exhibits  the  same  overall  variation  with  v  for  all  three 
boundaries.  If  the  flow  is  laminar,  (Ac)^  increases  linearly  with  mean  speed 
for  any  L,  whereas  for  the  mixed  flow,  the  increase  is  nonlinear.  This  can 
be  understood  from  the  following  considerations.  To  the  first-order  in  M,  a 
Taylor  expansion  of  ^ K\  —  DL  yields: 


(A c)L  ~ 


Dl 

A ! 


1  Cl 

A£  ci0)  ’ 


(39) 


where  Eq.  (31)  was  used.  Moreover,  the  norm  of  (Ci/C^^j  reflects  the  norm 
of  V,  and  for  the  laminar  flow,  this  norm  varies  directly  with  v.  Hence,  Eq.  (39) 
results  in  linear  variations  of  (A c)l  with  v  when  the  flow  is  laminar.  For  the 
mixed  flow,  however,  the  velocity  distribution  is  nonlinear  in  v  via  the  factor 
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a  (see  Eq.  (37)),  and  a  depends  nonlinearly  on  v  through  its  dependence  on 
Reynold’s  number  (cf.  Eq.  (38)).  The  nonlinearity  is  clearly  seen  in  all  the 
plots  for  the  mixed  flow  for  all  the  boundary  conditions  considered. 

The  second  distinguishing  feature  is  an  overall  fanning  out  of  (A c)l  as  s 
increases,  all  other  parameters  (£,  mean  flow,  frequency  and  the  waveguide) 
remaining  constant.  This  follows  from  the  fact  that  the  eigenvalues,  Kl ,  are 
monotonically  increasing  functions  of  s,  resulting  in  the  decrease  in  the  phase 
speed,  Ai.  In  addition,  the  spread  undergoes  a  general  increase  with  both  mode 
number,  £,  and  mean  flow.  The  increase  with  £  follows  from  the  interlacing 
property  [22]  of  the  zeros  of  Bessel  functions,  namely,  Kgti  <  Ke+i,i  <  < 

K£+it2  <  3  <  •  •  • .  Furthermore,  since  (C^ / C'|0)  j  ,  increases  with  v,  the  spread 

of  the  phase  speed  deviations  increases  with  the  mean  flow  speed. 

In  the  Robin  boundary,  attenuation  of  the  waves  may  be  expected.  In  other 
words,  (A c)l  may  have  imaginary  components.  It  can  be  seen  from  Eq.  (35) 
that  attenuation  is  expected  to  be  significant  near,  and  of  course,  above  the 
cut-off  frequency  of  a  mode.  For  the  modes  reported  in  this  paper,  the  numer¬ 
ical  calculations  showed  that  the  modes  were  away  from  their  cut-off  limits 
at  the  carrier  frequency  of  1  MHz.  Attenuation  was,  therefore,  not  expected 
to  be  significant.  The  results  of  the  numerical  computations  corroborate  this, 
as  can  be  seen  from  Fig.  4,  where  some  typical  variations  of  the  attenuation 
with  flow  are  plotted.  Figure  4(a)  show  Im(Ac)z,  vs.  v  for  (  =  0.3,  l  =  0, 
and  for  the  mixed  flow.  The  corresponding  results  for  (  =  1.2  are  plotted 
in  Fig.  4(b).  The  attenuation  is  seen  to  be  two  orders  of  magnitude  smaller 
than  Re(Ac)i.  Qualitatively  similar  behavior  was  observed  for  the  other  mode 
numbers.  Im(Ac)x,  for  l  —  2  in  Fig.  3(b)  are  smaller  in  magnitude  than  those 
for  l  =  0  in  Fig.  3(a).  This  behavior  was  typical  of  Im(Ac)^  irrespective  of 
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whether  the  flow  was  laminar  or  mixed. 


We  conclude  our  main  discussions  with  the  following  remark. 


Remark:  The  integral  Eq.  (8)  is  quite  general,  and  can  accommodate  practical 
physical  situations.  The  incident  field  that  appears  in  this  equation  can  be 
arbitrary,  and  the  medium  (or  a  part  thereof)  can  be  inhomogeneous.  In  the 
latter  case,  the  differential  Eq.  (4)  must  be  modified  to: 


2  i&Q. 


V2  +  /cq(1  +  e(r))  H - — V (r)  ■  V 


cu 


p(r)  =  0, 


(40) 


and  the  corresponding  inhomogeneous  part,  Pih(r),  becomes: 


Pih(r)  =  j  dr'G(0)(r|r) 
R 1 


o  7-p 

kle(j')X(v')  +  -^V(r')  •  V 
cu 


p(r). 


(41) 


In  Eqs.  (40)  and  (41),  e  is  the  medium’s  inhomogeneity,  and  x  is  the  charac¬ 
teristic  function  of  the  compact  region  occupied  by  e.  The  acoustic  velocity 
inhomogeneity  is  now  augmented  by  the  additional  inhomogeneity,  e(r'),  which 
provides  physical  scattering  centers.  Equations  (40)  and  (41)  are  representa¬ 
tives  of  numerous  practical  engineering  problems  involving  flow  mixed  with 
particulate  matters.  The  above  equations  hold  for  M  <  1.  In  addition,  if  e  is 
smaller  than  k0,  then  Born  approximation  can  be  invoked.  In  the  general  case, 
however,  e  may  not  be  small,  and  the  second  order  term  in  the  velocity  may 
have  to  be  incorporated.  Fortunately,  the  integral  in  Eq.  (41)  is  a  convolution 
integral,  and  can  be  accomplished  by  fast  solvers  such  as  CG-FFT  (conjugate 
gradient-fast  Fourier  transform)  [25,26].  Furthermore,  depending  upon  the 
relative  magnitudes  of  the  velocity  and  the  physical  scattering  centers,  WKB 
type  approximation  can  be  applied. 
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Fig.  3.  Re(Ac)^  vs.  v  for  the  Robin  boundary  conditions,  /  =  1MHz:  (a)  (  =  0.3, 
laminar  flow,  i  =  0,  (b)  £  =  0.3,  mixed  flow,  i  =  2  (c)  C  =  1-2,  laminar  flow,  i  =  0, 
(d)  C  =  1.2,  mixed  flow,  i  =  2.  solid:  s  =  0,  dot:  s  =  1,  dash:  s  =  2,  dash-dot:  s  =  3, 
dash-dot-dot-dot:  s  =  4 


((Ac),)  (m 


Fig.  4.  Irn(Ac),  vs.  v  for  mixed  flow  with  the  Robin  boundary  conditions,  £  =  1.2, 
/  =  1MHz:  (a)  l  =  1,  (b)  £  =  2.  dot:  s  =  1,  dash:  s  =  2,  dash-dot:  s  =  3, 
dash-dot-dot-dot:  s  =  4 
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